home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio / DSeqFile.cpp < prev    next >
Text File  |  1996-07-05  |  12KB  |  488 lines
  1. // DSeqFile.cp 
  2. // d.g.gilbert
  3.     
  4.  
  5. #include <ncbi.h>
  6. #include <dgg.h>
  7. #include <stdio.h>
  8. #include <DFile.h>
  9. #include "DSeqFile.h"
  10. #include "DSeqList.h"
  11. #include "ureadseq.h"
  12. #include "DMegaSequence.h"
  13. #include "DSequence.h"
  14.  
  15.  
  16.  
  17. Global char* gOutputName = NULL;
  18.  
  19. Boolean DSeqFile::gWriteMasks = true;
  20. Boolean DSeqFile::gSaveMasks= true;
  21. static char* kMaskName = "#Mask";
  22.  
  23. // static
  24. void DSeqFile::DontSaveMasks()
  25. {
  26.     gSaveMasks= gWriteMasks;
  27.     gWriteMasks= false;
  28. }
  29.  
  30. void DSeqFile::DoSaveMasks()
  31. {
  32.     gWriteMasks= gSaveMasks;
  33. }
  34.  
  35.  
  36. Boolean dupBases = false; // !?
  37.  
  38. //static
  39. DSequence* DSeqFile::MakeSequence(char* name, char*& bases, char* info, long modtime)
  40. {
  41.     char *b;
  42.     long nbases;
  43.     Boolean ismega;
  44.     DSequence* aSeq;
  45.     
  46.     for (b= bases, nbases=0; 
  47.             (b) && (*b) && nbases<DMegaSequence::kMegaPartSize; 
  48.             b++, nbases++)
  49.                 ;
  50.     ismega= nbases>=DMegaSequence::kMegaPartSize; 
  51.      
  52.     if (ismega) {
  53.         aSeq= new DMegaSequence();
  54.         } 
  55.     else {
  56.         aSeq= new DSequence();
  57.         }
  58.     aSeq->SetName(name);
  59.     aSeq->SetBases(bases, dupBases);
  60.     aSeq->SetInfo(info);
  61.      if (modtime!=0) aSeq->SetTime( (unsigned long)modtime);
  62.     aSeq->UpdateFlds();      
  63.     return aSeq;
  64. }
  65.  
  66.  
  67.  
  68. #if NEED_LATER
  69. // must be in scope of a DSeqFile:: method
  70. Local const  struct formatTable {
  71.   char  *name;
  72.   short num;
  73.   } formname[] = {
  74.     {"ig",  kIG},
  75.     {"stanford", kIG},
  76.     {"genbank", kGenBank},
  77.     {"gb", kGenBank},
  78.     {"nbrf", kNBRF},
  79.     {"embl", kEMBL},
  80.     {"gcg", kGCG},
  81.     {"uwgcg", kGCG},
  82.     {"dnastrider", kStrider},
  83.     {"strider", kStrider},
  84.     {"fitch", kFitch},
  85.     {"pearson", kPearson},
  86.     {"fasta", kPearson},
  87.     {"zuker", kZuker},
  88.     {"olsen", kOlsen},
  89.     {"phylip", kPhylip},
  90.     {"phylip3.2", kPhylip2},
  91.     {"phylip3.3", kPhylip3},
  92.     {"phylip3.4", kPhylip4},
  93.     {"phylip-interleaved", kPhylip4},
  94.     {"phylip-sequential", kPhylip2},
  95.     {"plain", kPlain},
  96.     {"raw", kPlain},
  97.     {"pir", kPIR},
  98.     {"codata", kPIR},
  99.     {"asn.1", kASN1},
  100.     {"msf", kMSF},
  101.     {"paup", kPAUP},
  102.     {"nexus", kPAUP},
  103.     {"pretty", kPretty},
  104.     {0,0}
  105.   };
  106. #endif
  107.  
  108. enum { kMaxFormatNames = 18 }; // != kMaxFormat ??
  109. static char *gSeqFileFormats[] = { //DSeqFile::kMaxFormat+1
  110.         "IG|Stanford",
  111.         "GenBank|GB",
  112.         "NBRF",
  113.         "EMBL",
  114.         "GCG",
  115.         "DNAStrider",
  116.         "Fitch",
  117.         "Pearson|Fasta",
  118.         "Zuker [in only]",
  119.         "Olsen [in only]",
  120.         "Phylip3.2|Phylip2",
  121.         "Phylip|Phylip4",
  122.         "Plain|Raw",
  123.         "PIR|CODATA",
  124.         "MSF",
  125.         "ASN.1",
  126.         "PAUP|NEXUS",
  127.         "Pretty [out only]",// # 18
  128.         //"AutoSeq [in only]", 
  129.         0};
  130.  
  131. const char* DSeqFile::FormatName(short index)
  132. { 
  133.     if (index<1 || index>kMaxFormatNames) return "";
  134.     else return gSeqFileFormats[index-1];
  135. }
  136.  
  137. short  DSeqFile::FormatFromName(const char* name)
  138. {
  139.     if (!name) return kUnknown;
  140.     long   i, len1, namelen= StrLen(name);
  141.     char      *form1, *form2;
  142.     if (namelen<2) return kUnknown;
  143.     for (i=0; i<kMaxFormat; i++) {
  144.         form1= (char*) gSeqFileFormats[i];
  145.         form2= StrChr(form1,'|'); 
  146.         if (form2) { len1= form2 - form1; form2++; }
  147.         else len1= StrLen( form1);
  148.         if (StrNICmp(name, form1,len1) == 0) return i+1;
  149.         else if (form2 && StrICmp(name, form2) == 0) return i+1;
  150.         }
  151.     return kUnknown;
  152. }
  153.  
  154.  
  155. char* DSeqFile::FixSeqID( char* s, short leftmargin)
  156. {
  157.     register char * cp, * ep;
  158.     for (cp= s; *cp && isspace(*cp); cp++) ;
  159.     for (ep= cp; *ep && !isspace(*ep); ep++) ;
  160.     if (ep>cp && ep[-1] == ',') ep--; // fix for fasta "name, #bases, # checksum"
  161.     if (*ep) *ep= 0;
  162.     return cp;
  163. }
  164.  
  165. short DSeqFile::SeqFileFormatWrapper( DFile* aFile) 
  166. {
  167.     long    skiplines = 0;
  168.     short    aFormat, err = 0;
  169.  
  170.     aFile->Open("r");  
  171.     aFormat= seqFileFormatFp( aFile->fFile, &skiplines, &err);
  172.     aFile->Close();  
  173.     return    aFormat;
  174. }
  175.  
  176.  
  177. short DSeqFile::ReadSeqFile( DFile* aFile, DSeqList* aSeqList) 
  178. {
  179.     char    seqid[256], *seqidptr, *bases, *seqlist;
  180.     long    skiplines = 0, skipi, seqlen;
  181.     short aFormat, atseq = 0, whichSeq = 0, nseq = 0, err = 0;
  182.     char  *defname;
  183.     Boolean notdone = true, needrewind = false;
  184.     DSequence* aSeq;
  185.     
  186.     aFile->OpenFile("r");    
  187.     aFormat = seqFileFormatFp( aFile->fFile, &skiplines, &err);
  188.     aFile->OpenFile(); // rewind
  189.     defname= FixSeqID( gDefSeqName, 0);
  190.     
  191.     switch (aFormat) {
  192.         case kUnknown:
  193.         case kNoformat:
  194.             aFile->CloseFile();
  195.             return aFormat;
  196.          
  197.             // all multi-pass (interleaved) formats go here !        
  198.         case kOlsen:
  199.         case kMSF:
  200.         case kPAUP:
  201.         case kPhylip2:
  202.         case kPhylip4:
  203.             needrewind = true;
  204.             aFile->CloseFile();
  205.             seqlist = listSeqs( aFile->GetName(), skiplines, aFormat, &nseq, &err);
  206.             MemFree(seqlist); 
  207.             break;
  208.         }
  209.         
  210.     skipi= skiplines;
  211.     while (notdone) {
  212.             // this one-pass thru gFileRef may well be buggy for current readSeq
  213.             // need to rewind file & count whichSeq++ for each sequence...?
  214.         if (needrewind) {
  215.             whichSeq++; 
  216.             atseq= 0; 
  217.             aFile->OpenFile("r");
  218.             skipi= skiplines;
  219.             }
  220.         else {
  221.             whichSeq = 1; 
  222.             atseq= 0;
  223.             }
  224.         seqlen = 0; seqidptr= seqid; *seqidptr= 0;
  225.         bases = readSeqFp( whichSeq, aFile->fFile, skipi, aFormat, &seqlen, 
  226.                                                         &atseq, &err, seqidptr);
  227.         skipi= 0; 
  228.         
  229.                 // install bases into new TSequence, and stick TSequence into list
  230.         if (bases) {
  231.             if (seqlen > 0) {
  232.                 char *name, * cp;
  233.                 Boolean ismask= false;
  234.                 
  235.                 if (*seqidptr != 0) {
  236.                     if ( (cp= StrStr(seqidptr, kMaskName)) != NULL) {
  237.                         ismask= true;
  238.                         *cp= 0;
  239.                         }
  240.                     name= FixSeqID( seqidptr, 0);
  241.                     }
  242.                 else 
  243.                     name= defname;
  244.                     
  245.                 if (ismask) {
  246.                     aSeq= aSeqList->FindNamedSeq( name);
  247.                     if (aSeq) {
  248.                         // convert mask...
  249.                         register char *mp;
  250.                         for (mp= bases; *mp; mp++) { 
  251.                             *mp -= '?'; // to zero relative
  252.                             *mp |= 0x80; // add high bit
  253.                             }
  254.                         aSeq->SetMasks( bases);
  255.                         }
  256.                     }
  257.                 else {
  258.                     aSeq = MakeSequence( name, bases, NULL, 0); //gFileRef->moddate
  259.                     aSeqList->InsertLast(aSeq);
  260.                     }
  261.                 }
  262.             else notdone= false;
  263.             MemFree(bases); // readSeq made this ptr w/ malloc call... MakeSeq copies this data... 
  264.             }
  265.         else notdone= false; // nodata -- if format=0, no data is read
  266.             
  267.         if (needrewind) notdone &= (whichSeq < nseq); // stop now
  268.         else notdone &= !aFile->Eof();  
  269.         }
  270.         
  271.     aFile->CloseFile();
  272.     return aFormat;
  273. }
  274.  
  275.  
  276.  
  277.  
  278. void    DSeqFile::WriteSeqHeader( DFile* aFile,  
  279.                                 short format, short nseqs, long nbases, char* outputName,
  280.                                 Boolean& needSameSize, Boolean& isInterleaved)
  281. {
  282.     const char *kASN1headline = "Bioseq-set ::= {\nseq-set {\n";
  283.     char    line[256], name[64];
  284.     
  285.     // this is a ureadseq.h call
  286.     // need to do somewhere once, then let user set prefs...
  287.  
  288. #define gPrettyInit1(p) { \
  289.   p.isactive=false;\
  290.   p.baseonlynum=true;\
  291.   p.numline= p.atseq= 0;\
  292.   p.numright= p.numtop= true; \
  293.   p.numleft= p.numbot= false;\
  294.   p.nameright= true; \
  295.   p.nameleft= p.nametop= false;\
  296.   p.noleaves= p.domatch= p.degap= false;\
  297.   p.matchchar='.';\
  298.   p.gapchar='-';\
  299.   p.namewidth=8;\
  300.   p.numwidth=5;\
  301.   p.interline=1;\
  302.   p.spacer=10;\
  303.   p.seqwidth=50;\
  304.   p.tab=0; }
  305.  
  306.     gPrettyInit1(gPretty);
  307.     //gPretty.interline= 2; // not active, along w/ numtop,numbot,nametop...
  308.     // need to work into this class
  309.     
  310.     
  311.         
  312.     needSameSize= false;
  313.     isInterleaved= false;
  314.     aFile->Open("a"); 
  315.  
  316.     if (outputName) StrNCpy( name, outputName, 63);
  317.     else StrCpy( name, DFileManager::kUntitled);
  318.     
  319.     switch (format) {
  320.         case kPIR:
  321.             aFile->WriteLine( "\\\\\\\n");
  322.             break;
  323.         case kASN1:
  324.             aFile->WriteLine( (char*)kASN1headline);
  325.             break;
  326.  
  327.         case kMSF:     // this is for top of interleaved file
  328.             {
  329.             isInterleaved= true;
  330.          unsigned long checksum = 0, checkall = 0;
  331.             //checksum= GCGchecksum(seq, seqlen, &checkall);
  332.             sprintf(line,"\n %s  MSF: %d  Type: N  June 01, 1776 12:00  Check: %d ..\n\n",
  333.                                         name, nbases, checkall);
  334.             aFile->WriteLine( line);
  335.             }
  336.             break;
  337.  
  338.         case kPhylip2:     
  339.             isInterleaved= false;
  340.             goto allphyl;
  341.             break;
  342.  
  343.         case kPhylip4:      
  344.             isInterleaved= true;
  345.         allphyl:
  346.             {
  347. //    if (phylvers >= 4) fprintf(foo," %d %d\n", seqout, seqlen);
  348. //    else fprintf(foo," %d %d YF\n", seqout, seqlen);
  349.             if (isInterleaved) sprintf( line, " %d %d\n", nseqs, nbases);  // " %d %d F\n"
  350.             else sprintf( line, " %d %d I \n", nseqs, nbases);  // " %d %d FI \n"
  351.             aFile->WriteLine( line);
  352.             needSameSize= true;
  353.             }
  354.             break;
  355.  
  356.         case kPAUP:
  357.             isInterleaved= true;
  358.             needSameSize= true;
  359.             aFile->WriteLine( "#NEXUS\n");
  360.             sprintf( line,"[%s -- data title]\n\n", name);
  361.             aFile->WriteLine( line);
  362.             break;
  363.  
  364.         case kPretty:
  365.             isInterleaved= true;
  366.             gPretty.isactive= true;
  367.             break;
  368.             
  369.         default:
  370.             break;
  371.         }
  372.         
  373.     aFile->CloseFile();
  374. }
  375.  
  376.  
  377.  
  378. void    DSeqFile::DoInterleave( DFile* toFile, DFile* fromFile, 
  379.                                 short outform, short seqtype, short nseqs, long nbases,
  380.                                 short nlines, short noutindex, long* outindex)
  381. {
  382.     short        leaf, iline, iseq;
  383.     char        *cp, line[256];
  384.     
  385.     toFile->Open("a");
  386.     fromFile->Open("r");
  387.     
  388.     for (leaf=0; leaf<nlines; leaf++) {
  389.     
  390.         if (outform == kMSF && leaf == 1) { // do after writing seq names
  391.             toFile->WriteLine( "//\n\n");
  392.             }
  393.             
  394.         if (outform == kPAUP && leaf==1) { // do after writing seq names
  395.             switch (seqtype) {
  396.                 case DSequence::kDNA     : cp= "dna"; break;
  397.                 case DSequence::kRNA     : cp= "rna"; break;
  398.                 case DSequence::kNucleic : cp= "dna"; break;
  399.                 case DSequence::kAmino   : cp= "protein"; break;
  400.                 case DSequence::kOtherSeq: cp= "dna"; break;
  401.                 }
  402.             toFile->WriteLine( "\nbegin data;\n");
  403.             sprintf(line," dimensions ntax=%d nchar=%d;\n", nseqs, nbases); 
  404.             toFile->WriteLine( line);
  405.             sprintf(line," format datatype=%s interleave missing=%c", cp, gPretty.gapchar); // gapchar == indelHard ??? 
  406.             toFile->WriteLine( line);
  407.             //if (gPretty.domatch) { sprintf(line," matchchar=%c", gPretty.matchchar); 
  408.                 //Mwrite( toFile, line); }
  409.             toFile->WriteLine( ";\n  matrix\n");
  410.             }
  411.  
  412.         for (iseq=0; iseq<noutindex; iseq++) {
  413.             fromFile->SetDataMark( outindex[iseq]);
  414.             
  415.             for (iline=0; iline<=leaf; iline++)
  416.                 if (fromFile->ReadLine( line, 256) != 0) *line= 0;
  417.                 
  418.             if (fromFile->Tell() <= outindex[iseq+1])
  419.                 toFile->WriteLine( line, false);  //remember- ReadLine a.k.a fgets retains newline
  420.             }
  421.  
  422.      for (iline=0; iline<gPretty.interline; iline++)
  423.         toFile->WriteLine("\n"); //LineEnd 
  424.             
  425.         }
  426.         
  427.     fromFile->Close();
  428.     toFile->Close();
  429. }
  430.  
  431.  
  432. void    DSeqFile::WriteSeqTrailer( DFile* aFile,
  433.                                 short format, short /*nseqs*/, short /*nlinesperseq*/, long /*nbases*/)
  434. {
  435.     switch (format) {
  436.   
  437.       case kASN1:  
  438.             aFile->Open("a"); 
  439.             aFile->WriteLine("} }\n"); 
  440.             aFile->Close(); 
  441.             break;
  442.     
  443.       case kPAUP:    
  444.             aFile->Open("a"); 
  445.             aFile->WriteLine(";\n  end;\n");
  446.             aFile->Close(); 
  447.             break;    
  448.         }
  449. }
  450.  
  451.  
  452.  
  453. void DSeqFile::WriteSeqWrapper( DFile* aFile,
  454.                                     char *seq, long seqlen, short outform, char *seqid)
  455. { 
  456.     if (!aFile->IsOpen())  aFile->Open("a"); 
  457.     seqid= StrDup(seqid);
  458.     char* cp= FixSeqID(seqid, 0);
  459.     gLinesPerSeqWritten= writeSeq( aFile->fFile, seq, seqlen, outform, cp);
  460.     MemFree(seqid);
  461.     //aFile->Close(); //?? want left open for appending by caller ? e.g., seqmail
  462. }  
  463.  
  464. void DSeqFile::WriteMaskWrapper( DFile* aFile,
  465.                                     char *mask, long seqlen, short outform, char *seqid)
  466. { 
  467.     char maskid[128];
  468.     
  469.     if (!aFile->IsOpen())  aFile->Open("a"); 
  470.     seqid= StrDup(seqid);
  471.     char* cp= FixSeqID(seqid, 0);
  472.     StrNCpy( maskid, cp, sizeof(maskid)-5);
  473.     StrNCat( maskid, kMaskName, sizeof(maskid));
  474.     MemFree(seqid);
  475.         //     massage mask into writable form
  476.     mask= StrDup(mask); 
  477.     Boolean hasdata= false;
  478.     for (cp= mask; *cp; cp++) { 
  479.         *cp &= 0x7F; // drop hi bit
  480.         if (*cp != 0) hasdata= true;
  481.         *cp += '?';  // set 0 == '?', mask range is mostly 'A-Za-z'
  482.      }
  483.     if (hasdata)
  484.         gLinesPerSeqWritten= writeSeq( aFile->fFile, mask, seqlen, outform, maskid);
  485.     MemFree(mask);
  486. }  
  487.  
  488.